# Automatic network construction, module detection and module relating trains.

library(WGCNA)

args <- commandArgs(trailingOnly = TRUE)

exprFilePath = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/datExpr.csv"
traitFilePath = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/datTraits.csv"

datExpr = read.csv(exprFilePath, stringsAsFactors = FALSE, row.names = 1);
datTraits = read.csv(traitFilePath, stringsAsFactors = FALSE, row.names = 1);

nGenes = ncol(datExpr);
nSamples = nrow(datExpr);

# construct network
net = blockwiseModules(datExpr, power = 6,
                       TOMType = "unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/TOM",
                       verbose = 3)

# save net data
save(net, file="network.RData")

# Visualizatioin network modules
png(
  filename = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/network_module.png",
  width = 13,
  height = 9,
  units = "in",
  res = 300
)
par(mfrow = c(1,2));
cex1 = 0.9;
mergedColors = labels2colors(net$colors)
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

dev.off()


# Visualizatioin network heatmap
png(
  filename = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/network_heatmap.png",
  width = 13,
  height = 9,
  units = "in",
  res = 300
)

moduleColors = labels2colors(net$colors)
dissTOM = 1-TOMsimilarityFromExpr(datExpr, power = 6);

# 随机选择部分基因
nSelect = ceiling(nGenes / 10)
set.seed(10);
select = sample(nGenes, size = nSelect);
selectTOM = dissTOM[select, select];
selectTree = hclust(as.dist(selectTOM), method = "average")
selectColors = moduleColors[select];
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes")

dev.off()


# Relate modules to datTraits
nGenes = ncol(datExpr);
nSamples = nrow(datExpr);
# Recalculate MEs with color labels
moduleColors = labels2colors(net$colors)
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
MEs = orderMEs(MEs0)
moduleTraitCor = cor(MEs, datTraits, use = "p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples);
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
                   signif(moduleTraitPvalue, 1), ")", sep = "");

png(
  filename = "/Users/zyl/Workspace/hzau/MOIN-Extension-Github/WGCNA-Extension/FemaleLiver-Data/module_trait.png",
  width = 13,
  height = 9,
  units = "in",
  res = 300
)
par(mar = c(6, 8.5, 3, 3));
labeledHeatmap(Matrix = moduleTraitCor,
               xLabels = names(datTraits),
               yLabels = names(MEs),
               ySymbols = names(MEs),
               colorLabels = FALSE,
               colors = blueWhiteRed(50),
               textMatrix = textMatrix,
               setStdMargins = FALSE,
               cex.text = 0.5,
               zlim = c(-1,1),
               main = paste("Module-trait relationships"))
dev.off()
